home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 212_01 / percolar.c < prev    next >
Encoding:
C/C++ Source or Header  |  1980-01-01  |  18.3 KB  |  783 lines
  1. /* PERCOLAR.C    VERS:- 01.00  DATE:- 10/06/86  TIME:- 05:23:14 PM */
  2. /*
  3. %CC1 $1.C -X -O -E 7000
  4. %CLINK $1 DIO -F FLOAT -S
  5. %DELETE $1.CRL 
  6. */
  7. /* 
  8. Description:
  9.  
  10. Simulation of percolation on a two-dimensional square lattice, 
  11. Option to block one or more regions of the lattice to prevent sites 
  12. being filled (approximately circular blocked region = 1,3,5,3,1).
  13. The algorithm is not efficient (300 seconds/1000 site 2-d sheet).
  14.  
  15. A percolative process is characterized by a percolation transition, a phase
  16. change at which long-range connectivity (e.g., a gel phase; conduction 
  17. across the sample; etc) suddenly appears.  Connectivity is defined by the
  18. extent of a cluster of filled sites adjacent (ie, connected) in the lattice.
  19. Long-range connectivity is equivalent to a cluster of infinite extent.
  20.  
  21. The percolation model applies to stochastic processes taking place 
  22. in a spatially disordered system.  
  23. Examples: 
  24.     gelation; 
  25.     flow in a porous medium; 
  26.     dilute magnets; 
  27.     conductivity of conductor-insulator composite materials.
  28.  
  29. For a discussion of percolation, see R. Zallen, "The Physics of Amorphous
  30. Solids", John Wiley, New York, 1983.
  31.  
  32. This program finds the fractional occupancy of the model lattice at the 
  33. percolation transition.  Lattice sites are filled randomly.  When each filled
  34. site is added to the lattice, connectivity is tested along one dimension
  35. of the lattice (ie, between one pair of arbitrarily chosen sides of the 
  36. lattice).  The point of the percolation transition is the point at which
  37. connectivity between sides first appears.
  38.  
  39. Because of the stochastic character of percolation, a simulation consists
  40. of the averaging of NRPT determinations of the percolation transition.
  41.  
  42. The program is set for looping over the simulation, with change in
  43. control parameters at each loop.
  44.  
  45. Requires sqrt() function.
  46.  
  47. By J.A. Rupley, Tucson, Arizona
  48. Coded for BDS C compiler, version 1.50a
  49. */
  50.  
  51. #include "jario.h"
  52. #include "bdscio.h"
  53. #include "dio.h"
  54.  
  55. /*CHANGE TO FALSE IF SYSTEM NOT CONFIGURED FOR HAYES MODEM AS PUNCH-READER*/
  56. #define HAYES_MODEM FALSE
  57.  
  58. #DEFINE ARRAY_DIM 75
  59. #DEFINE NROW  10
  60. #DEFINE NCOL  10
  61. #DEFINE NRPT  20
  62. #DEFINE NBLOCKED 0
  63. #DEFINE NO_PRINT 2
  64.  
  65. #DEFINE PTS_BLOCKED 13        /*reset if change blocking pattern*/
  66.  
  67. #DEFINE FILLED '+'
  68. #DEFINE BLOCKED 'B'
  69. #DEFINE TRAIL 'O'
  70. #DEFINE BLANK ' '
  71.  
  72. #DEFINE ROW_PATTERN 2;
  73. #DEFINE COL_PATTERN 7;
  74. #DEFINE CIRCLE TRUE;
  75. #DEFINE BAR !CIRCLE;
  76.  
  77. int sum[ARRAY_DIM];
  78. char mast_array[ARRAY_DIM][ARRAY_DIM];
  79. char work_array[ARRAY_DIM][ARRAY_DIM];
  80.  
  81. int cnt_crit[100], cnt_blocked[100];
  82.  
  83. char del_i[8], del_j[8];
  84. char test_block[30];
  85. char str_buf[80];
  86. char title[160];
  87.  
  88. int xnrpt, xnrow, xncol;
  89. int nrpt, nrpt_1;
  90. int nsites, nsites_blocked;
  91. int no_print;
  92. int nblocked;
  93. int circle, bar, row_pattern, col_pattern;
  94.  
  95. char i_last, j_last, k_last;
  96.  
  97. int block_sum, crit_sum;
  98. char fcrit_avg[5], fcrit_var[5], fcrit_1var[5];
  99. char fract_avg[5], fract_var[5], fract_1var[5];
  100. char fbloc_avg[5], fbloc_var[5], fbloc_1var[5];
  101. char frat1_avg[5], frat1_var[5], frat1_1var[5];
  102. char frat2_avg[5], frat2_var[5], frat2_1var[5];
  103. char funblocked[5], fnrpt[5], fnrpt_1[5], ftemp[5], fsites[5], fsqd[5];
  104.  
  105. char i, j, k, ii, jj, kk, ll;
  106. int count;
  107. char pspecial;
  108. int s1, s2, s3;
  109. char clock_in[12];
  110.  
  111. /* page eject*/
  112.  
  113.  
  114. main(argc, argv)
  115. int argc;
  116. char **argv;
  117. {
  118.     dioinit(&argc, argv);
  119.     /*PSPECIAL - LOOP OVER VARIED PARMS*/
  120.     /*SETUP*/
  121.     for (pspecial = 0; pspecial < 100; pspecial++)
  122.     {
  123.         switch (pspecial)
  124.         {
  125.         case 0 :        /*case 0 used for setup*/
  126.             setup();
  127.             /*insert title(up to 2 lines)*/
  128.             strcpy(title,
  129.             "SIMULATION OF PERCOLATION ON A TWO-DIMENSIONAL SQUARE LATTICE"
  130.                 );
  131.             strcat(title,
  132.             ""
  133.                 );
  134.             printf("\n%s\n", title);
  135.             read_variables();        /*********optionally delete*/
  136.             continue;
  137.         case 1 :        /*case 1 to 99 = specials*/
  138.             circle = FALSE;
  139.             bar = !circle;
  140.             break;
  141.         default :
  142.             dioflush();
  143.             exit();
  144.         }
  145.  
  146.         nsites = xnrow * xncol;
  147.         header_disp();
  148.  
  149.         /*MAIN LOOP*/
  150.         /*repeat determination of critical count */
  151.         /*to obtain the average*/
  152.         for (nrpt = 0; nrpt < xnrpt; nrpt++)
  153.         {
  154.             /*initialize mast_array */
  155.             /*with optional blocking of sites*/
  156.             init_array();
  157.  
  158.             /*SUB LOOP*/
  159.             /*start of cycle of filling sites,*/
  160.             /*until have connection top-to-bottom*/
  161.             /*the critical count is the number filled*/
  162.             /*at the point of first connection*/
  163.             while (TRUE)
  164.             {
  165.                 /*find unfilled site, by random search*/
  166.                 i = (nrand(1) % xnrow) + 1;
  167.                 j = (nrand(1) % xncol) + 1;
  168.                 if (mast_array[i][j] != BLANK)
  169.                     continue;
  170.                 /*add new site to the master array*/
  171.                 i_last = i;
  172.                 j_last = j;
  173.                 mast_array[i][j] = FILLED;
  174.                 if (exit_test())
  175.                     array_disp();
  176.                 /*if connection not possible,*/
  177.                 /*go fill another site*/
  178.                 if ((cnt_crit[nrpt] = sum_filled(mast_array)) == 0)
  179.                     continue;
  180.  
  181.                 /*test for connection*/
  182.                 /*if no connection, go fill another site*/
  183.                 copy_array(mast_array, work_array);
  184.                 if (!work_test())
  185.                     continue;
  186.  
  187.                 /*EXIT SUBLOOP*/
  188.                 break;
  189.             }
  190.  
  191.             /*display mast_array and work_array,*/
  192.             /*side by side*/
  193.             /*display blocked sites, if present*/
  194.             /*display cnt_crit values and their*/
  195.             /*average and variance*/
  196.             if (nrpt < no_print)
  197.             {
  198.                 array_disp();
  199.                 crit_calc();
  200.                 crit_disp();
  201.             }
  202.             else
  203.                 printf("working on cycle %d\n", (nrpt + 1));
  204.         }
  205.         /*ALL DONE--TERMINATION OF MAIN LOOP*/
  206.         prog_exit();
  207.     }
  208.     /*END OF PSPECIAL LOOP*/
  209. }
  210. /* END OF MAIN                */
  211. /* page eject*/
  212.  
  213.  
  214. int read_variables()
  215. {
  216.     printf
  217.         ("\nDo you want to change values of program variables? (y/cr=>default): ");
  218.     getline(str_buf, 2);
  219.     if ((str_buf[0] != 0x79) && (str_buf[0] != 0x59))
  220.     {
  221.         printf("\n");
  222.         return;
  223.     }
  224.     printf
  225.         ("\nNew title? \n*:");
  226.     getline(str_buf, 79);
  227.     if (strlen(str_buf))
  228.     {
  229.         strcpy(title, str_buf);
  230.         printf("     Another title line?\n*:");
  231.         getline(str_buf, 79);
  232.         if (strlen(str_buf))
  233.             strcat(title, "\n");
  234.         strcat(title, str_buf);
  235.     }
  236.     printf("\n");
  237.     read_val(
  238.     "How many array rows? (value/cr=>no change = %d): ",
  239.     &xnrow);
  240.     printf("\n");
  241.     read_val(
  242.     "How many array cols? (value/cr=>no change = %d): ",
  243.     &xncol);
  244.     printf("\n");
  245.     read_val(
  246.     "How many blocked regions? (value/cr=>no change = %d): ",
  247.     &nblocked);
  248.     printf("\n");
  249.     read_val(
  250.     "How many repeat determinations? (value/cr=>no change = %d): ",
  251.     &xnrpt);
  252.     printf("\n");
  253.     read_val(
  254.     "How many cycles do you want printed? (value/cr=>no change = %d): ",
  255.     &no_print);
  256.  
  257.     printf("\n\n");
  258.     return;
  259. }
  260. /*END OF READ_VARIABLES            */
  261.  
  262.  
  263.  
  264. int read_val(string, val)
  265. char *string;
  266. int *val;
  267. {
  268.     while (TRUE)
  269.     {
  270.         printf(string, *val);
  271.         if (strlen(gets(str_buf)) == 0)
  272.             return;
  273.         *val = atoi(str_buf);
  274.     }
  275.     return;
  276. }
  277. /*END OF READ_VAL            */
  278.  
  279.  
  280.  
  281. int setup()
  282. {
  283.     /* set default values of xnrow, etc*/
  284.     xnrow = NROW;
  285.     xncol = NCOL;
  286.     xnrpt = NRPT;
  287.     nblocked = NBLOCKED;
  288.     no_print = NO_PRINT;
  289.     nsites_blocked = PTS_BLOCKED;
  290.     row_pattern = ROW_PATTERN;
  291.     col_pattern = COL_PATTERN;
  292.     circle = CIRCLE;
  293.     bar = BAR;
  294.  
  295.     /*establish kernel for random number generator*/
  296.     read_clock();
  297.     nrand(-1, s1, s2, s3);
  298.  
  299.     /*fill arrays for choosing step (square lattice) in work_test()*/
  300.     initb(del_i, "-1,0,1,0,-1,0,1,0");
  301.     initb(del_j, "0,-1,0,1,0,-1,0,1");
  302.  
  303.     /*fill array used to set block of points*/
  304.     /*to be excluded from the fitting*/
  305.     initb(&test_block[0], "0,0,1,0,0");
  306.     initb(&test_block[5], "0,1,1,1,0");
  307.     initb(&test_block[10], "1,1,1,1,1");
  308.     initb(&test_block[15], "0,1,1,1,0");
  309.     initb(&test_block[20], "0,0,1,0,0");
  310.     for (i = 0; i < 25; i++)
  311.         if (test_block[i])
  312.             test_block[i] = BLOCKED;
  313.         else
  314.             test_block[i] = BLANK;
  315.  
  316.     return;
  317. }
  318. /*END OF SETUP                */
  319.  
  320.  
  321.  
  322. #if HAYES_MODEM
  323. int read_clock()
  324. /*set kernel for random number generator based on read of time/date/day*/
  325. {
  326.     /* Clock_in gives hhmmssx    */
  327.     hayes_chrono("RT\r", clock_in);
  328.     s1 = atoi(&clock_in[2]);
  329.     printf("\nTime, date and day =  %s ", clock_in);
  330.  
  331.     /* Clock_in gives yymmdd    */
  332.     hayes_chrono("RD\r", clock_in);
  333.     s2 = atoi(&clock_in[2]);
  334.     printf("%s ", clock_in);
  335.  
  336.     /* Clock_in gives  d = acsii #    */
  337.     hayes_chrono("RW\r", clock_in);
  338.     s3 = atoi(clock_in);
  339.     printf("%s \n", clock_in);
  340.  
  341.     /*set generator*/
  342.  
  343.     return;
  344. }
  345. /*END OF READ_CLOCK            */
  346.  
  347.  
  348.  
  349. #define PUNCH 6
  350. #define READER 7
  351. #define CR 0x0d
  352.  
  353. int hayes_chrono(command, buffer)
  354. /* Send command to clock and return with string*/
  355. char *command, *buffer;
  356. {
  357.     strcpy(str_buf, "AT");
  358.     strcat(str_buf, command);        /* Form command */
  359.  
  360.     for (i = 0; str_buf[i] != NULL; bios(PUNCH, str_buf[i++]))
  361.         ;
  362.  
  363.     for (i = 0; i < 12; i++)
  364.         if ((str_buf[i] = bios(READER, 0)) == CR)
  365.             break;
  366.     str_buf[i] = '\0';
  367.  
  368.     strcpy(buffer, str_buf);
  369.     return;
  370. }
  371. /* END OF HAYES_CHRONO    */
  372.  
  373.  
  374.  
  375. #else
  376. void read_clock()
  377. {
  378.     printf("\n\nPlease hit three keys, not too quickly, to seed random number generator\n\n");
  379.     s1 = s2 = s3 = 0;
  380.     while (!bios(2, 0))
  381.         s1++;
  382.     bios(3, 0);
  383.     while (!bios(2, 0))
  384.         s2++;
  385.     bios(3, 0);
  386.     while (!bios(2, 0))
  387.         s3++;
  388.     while (bios(2, 0))
  389.         bios(3, 0);
  390. }
  391. /* END OF "READ_CLOCK" FOR SYSTEM WITHOUT HAYES MODEM*/
  392. #endif
  393.  
  394.  
  395.  
  396. int init_array()
  397. {
  398.     /*loop until success or 100 times*/
  399.     ll = 100;
  400.     while (ll--)
  401.     {
  402.         /*zero arrays*/
  403.         for (i = 0; i < ARRAY_DIM; i++)
  404.         {
  405.             sum[i] = 0;
  406.             for (j = 0; j < ARRAY_DIM; j++)
  407.                 mast_array[i][j] = work_array[i][j] = BLANK;
  408.         }
  409.  
  410.         /*exclude block of sites*/
  411.         /*count excluded sites*/
  412.         cnt_crit[nrpt] = cnt_blocked[nrpt] = 0;
  413.         if (!nblocked)
  414.             return;
  415.         for (k = 0; k < nblocked; k++)
  416.         {
  417.             while (bar)
  418.             {
  419.                 i = (nrand(1) % xnrow) + 1;
  420.                 j = (nrand(1) % xncol) + 1;
  421.                 if ((i > 0) && (i < (xnrow - row_pattern + 2))
  422.                     && (j > 0)
  423.                     && (j < (xncol - col_pattern + 2)))
  424.                 {
  425.                     for (ii = (i); ii < (i + row_pattern); ii++)
  426.                         for (jj = (j); jj < (j + col_pattern); jj++)
  427.                         {
  428.                             if (mast_array[ii][jj] != BLOCKED)
  429.                             {
  430.                                 mast_array[ii][jj] = BLOCKED;
  431.                                 cnt_blocked[nrpt]++;
  432.                             }
  433.                         }
  434.                     break;
  435.                 }
  436.             }
  437.             while (circle)
  438.             {
  439.                 i = (nrand(1) % xnrow) + 1;
  440.                 j = (nrand(1) % xncol) + 1;
  441.                 if ((i > 2) && (i < xnrow - 1)
  442.                     && (j > 2) && (j < xncol - 1))
  443.                 {
  444.                     kk = 0;
  445.                     for (ii = (i - 2); ii < (i + 3); ii++)
  446.                         for (jj = (j - 2); jj < (j + 3); kk++, jj++)
  447.                             if ((mast_array[ii][jj] != BLOCKED)
  448.                                 && (test_block[kk] == BLOCKED))
  449.                             {
  450.                                 mast_array[ii][jj] = BLOCKED;
  451.                                 cnt_blocked[nrpt]++;
  452.                             }
  453.                     break;
  454.                 }
  455.             }
  456.         }
  457.         /*test for impermeable barrier to connection*/
  458.         for (i = 1; i < (xnrow + 1); i++)
  459.             for (j = 1; j < (xncol + 1); j++)
  460.                 if (mast_array[i][j] != BLOCKED)
  461.                     work_array[i][j] = FILLED;
  462.         if (work_test())
  463.             return;
  464.     }
  465.     printf("\nFailure in array initialization\n");
  466.     prog_exit();
  467.     dioflush();
  468.     exit();
  469. }
  470. /*END OF INIT_ARRAY            */
  471.  
  472.  
  473.  
  474. int exit_test()
  475. {
  476.     if (kbhit())
  477.     {
  478.         getline(str_buf, 2);
  479.         printf("\ntype x or X to exit or any other character to continue: ");
  480.         getline(str_buf, 2);
  481.         if ((str_buf[0] == 0x58) || (str_buf[0] == 0x78))
  482.         {
  483.             prog_exit();
  484.             dioflush();
  485.             exit();
  486.         }
  487.         return TRUE;
  488.     }
  489.     return FALSE;
  490. }
  491. /*END OF EXIT_TEST            */
  492.  
  493.  
  494.  
  495. int prog_exit()
  496. {
  497.     header_disp();
  498.     nrpt = nrpt - 1;
  499.     crit_calc();
  500.     crit_disp();
  501.     printf("\n\n\n\f");
  502.     return;
  503. }
  504. /*END OF PROG_EXIT            */
  505.  
  506.  
  507.  
  508. int header_disp()
  509. {
  510.     printf("\n\n%s\n", title);
  511. #if HAYES_MODEM
  512.     read_clock();
  513. #endif
  514.     printf
  515.         ("\nSimulation of percolation on a two-dimensional square lattice.");
  516.     printf
  517.         ("\nThe array size is %d rows by %d columns, %d sites total.",
  518.     xnrow, xncol, nsites);
  519.     printf
  520.         ("\nThere are %d regions of %d sites blocked from filling.\n",
  521.     nblocked, nsites_blocked);
  522.     return;
  523. }
  524. /*END OF HEADER_DISP            */
  525.  
  526.  
  527.  
  528. int array_disp()
  529. /*display master array and working array side-by-side, if possible*/
  530. /*display blocked sites, in mast_array*/
  531. /*display trail of the work_test pointer, in work_array*/
  532. /*display coordinates of last site filled*/
  533. {
  534.     printf("\n");
  535.     for (i = 1; i < (xnrow + 1); i++)
  536.     {
  537.         printf("\n");
  538.         for (j = 1; j < (xncol + 1); j++)
  539.             printf("%c", mast_array[i][j]);
  540.         if (xncol < 38)
  541.         {
  542.             printf("   ");
  543.             for (j = 1; j < (xncol + 1); j++)
  544.                 printf("%c", work_array[i][j]);
  545.         }
  546.     }
  547.     if (xncol > 37)
  548.     {
  549.         printf("\n\n\n");
  550.         for (i = 1; i < (xnrow + 1); i++)
  551.         {
  552.             printf("\n");
  553.             for (j = 1; j < (xncol + 1); j++)
  554.                 printf("%c", work_array[i][j]);
  555.         }
  556.     }
  557.     printf("\n\n nrpt=%d   cnt_crit[nrpt]=%d   ",
  558.     nrpt, cnt_crit[nrpt]);
  559.     printf("          i_last=%d   j_last=%d\n", i_last, j_last);
  560.     return;
  561. }
  562. /* END OF ARRAY_DISP            */
  563.  
  564.  
  565.  
  566. int crit_disp()
  567. /*display critical counts and fractions--averages and variances*/
  568. {
  569.     char fpno1[5], fpno2[5];
  570.     int sign;
  571.     char *sqrt();
  572.  
  573.     if (nblocked)
  574.     {
  575.         printf("\nList of all counts of blocked sites\n");
  576.         for (i = 0; i < (nrpt + 1); i++)
  577.             printf("%4d%c", cnt_blocked[i],
  578.             ((i % 10) == 9 || i == nrpt) ? '\n' : ' ');
  579.     }
  580.  
  581.     printf
  582.         ("\nList of all counts of filled sites at the critical point\n");
  583.     for (i = 0; i < (nrpt + 1); i++)
  584.         printf("%4d%c", cnt_crit[i],
  585.         ((i % 10) == 9 || i == nrpt) ? '\n' : ' ');
  586.  
  587.     printf(
  588.     "\n\n   critical count =%10.2f     sample_sd=%10.2f   avg_sd=%10.2f",
  589.     fcrit_avg,
  590.     sqrt(fpno1, &sign, fcrit_var),
  591.     sqrt(fpno2, &sign, fcrit_1var));
  592.     printf(
  593.     "\n   critical/total =%10.5f     sample_sd=%10.7f   avg_sd=%10.7f",
  594.     fract_avg,
  595.     sqrt(fpno1, &sign, fract_var),
  596.     sqrt(fpno2, &sign, fract_1var));
  597.  
  598.     if (nblocked)
  599.     {
  600.         printf(
  601.         "\n\n   blocked sites  =%10.2f     sample_sd=%10.2f   avg_sd=%10.2f",
  602.         fbloc_avg,
  603.         sqrt(fpno1, &sign, fbloc_var),
  604.         sqrt(fpno2, &sign, fbloc_1var));
  605.         printf(
  606.         "\n   blocked/total  =%10.5f     sample_sd=%10.7f   avg_sd=%10.7f",
  607.         frat1_avg,
  608.         sqrt(fpno1, &sign, frat1_var),
  609.         sqrt(fpno2, &sign, frat1_1var));
  610.         printf(
  611.         "\n   critical/unblkd=%10.5f     sample_sd=%10.7f   avg_sd=%10.7f",
  612.         frat2_avg,
  613.         sqrt(fpno1, &sign, frat2_var),
  614.         sqrt(fpno2, &sign, frat2_1var));
  615.     }
  616.     printf("\n\n\n");
  617.     return;
  618. }
  619. /* END OF CRIT_DISP            */
  620.  
  621.  
  622.  
  623. int crit_calc()
  624. /*calculate averages and variances of critical counts and fractions*/
  625. {
  626.     nrpt_1 = nrpt + 1;
  627.     itof(fnrpt, nrpt);
  628.     itof(fnrpt_1, nrpt_1);
  629.     itof(fsites, nsites);
  630.     fpmult(fsqd, fsites, fsites);
  631.  
  632.     crit_sum = 0;
  633.     for (i = 0; i < nrpt + 1; i++)
  634.         crit_sum = crit_sum + cnt_crit[i];
  635.  
  636.     itof(fcrit_avg, crit_sum);
  637.     fpdiv(fcrit_avg, fcrit_avg, fnrpt_1);
  638.     atof(fcrit_var, "0");
  639.     if (nrpt > 0)
  640.         for (i = 0; i < (nrpt + 1); i++)
  641.         {
  642.             itof(ftemp, cnt_crit[i]);
  643.             fpsub(ftemp, ftemp, fcrit_avg);
  644.             fpmult(ftemp, ftemp, ftemp);
  645.             fpadd(fcrit_var, ftemp, fcrit_var);
  646.         }
  647.     fpdiv(fcrit_var, fcrit_var, fnrpt);
  648.     fpdiv(fcrit_1var, fcrit_var, fnrpt_1);
  649.     fpdiv(fract_avg, fcrit_avg, fsites);
  650.     fpdiv(fract_var, fcrit_var, fsqd);
  651.     fpdiv(fract_1var, fcrit_1var, fsqd);
  652.  
  653.     /*calculations if blocked regions present*/
  654.     if (!nblocked)
  655.         return;
  656.     block_sum = 0;
  657.     for (i = 0; (i < (nrpt + 1)) && (cnt_blocked[i] != 0); i++)
  658.         block_sum = cnt_blocked[i] + block_sum;
  659.  
  660.     itof(fbloc_avg, block_sum);
  661.     fpdiv(fbloc_avg, fbloc_avg, fnrpt_1);
  662.     atof(fbloc_var, "0");
  663.     if (nrpt > 0)
  664.         for (i = 0; i < (nrpt + 1); i++)
  665.         {
  666.             itof(ftemp, cnt_blocked[i]);
  667.             fpsub(ftemp, ftemp, fbloc_avg);
  668.             fpmult(ftemp, ftemp, ftemp);
  669.             fpadd(fbloc_var, ftemp, fbloc_var);
  670.         }
  671.     fpdiv(fbloc_var, fbloc_var, fnrpt);
  672.     fpdiv(fbloc_1var, fbloc_var, fnrpt_1);
  673.  
  674.     fpdiv(frat1_avg, fbloc_avg, fsites);
  675.  
  676.     fpdiv(frat1_var, fbloc_var, fsqd);
  677.     fpdiv(frat1_1var, fbloc_1var, fsqd);
  678.  
  679.     fpsub(funblocked, fsites, fbloc_avg);
  680.     fpdiv(frat2_avg, fcrit_avg, funblocked);
  681.  
  682.     fpmult(frat2_var, frat2_avg, frat2_avg);
  683.     fpmult(frat2_var, frat2_var, fbloc_var);
  684.     fpadd(frat2_var, frat2_var, fcrit_var);
  685.     fpdiv(frat2_var, frat2_var, funblocked);
  686.     fpdiv(frat2_var, frat2_var, funblocked);
  687.  
  688.     fpmult(frat2_1var, frat2_avg, frat2_avg);
  689.     fpmult(frat2_1var, frat2_1var, fbloc_1var);
  690.     fpadd(frat2_1var, frat2_1var, fcrit_1var);
  691.     fpdiv(frat2_1var, frat2_1var, funblocked);
  692.     fpdiv(frat2_1var, frat2_1var, funblocked);
  693.  
  694.     return;
  695. }
  696. /*END OF CRIT_CALC            */
  697.  
  698.  
  699. int sum_filled(array)
  700. char array[ARRAY_DIM][ARRAY_DIM];
  701. /*sum filled sites of array, by row*/
  702. /*return 0 if any row is zero (which means there can be no connection
  703.     between top and bottom edges (rows))*/
  704. /*else return count of all filled sites in array*/
  705. {
  706.     count = 0;
  707.     for (i = 1; i < (xnrow + 1); i++)
  708.     {
  709.         sum[i] = 0;
  710.         for (j = 1; j < (xncol + 1); j++)
  711.             if (array[i][j] == FILLED)
  712.                 sum[i] = sum[i] + 1;
  713.         if (sum[i] == 0)
  714.             return 0;
  715.         count = count + sum[i];
  716.     }
  717.     return count;
  718. }
  719. /* END OF SUM_FILLED            */
  720.  
  721.  
  722.  
  723. int copy_array(source, dest)
  724. char source[ARRAY_DIM][ARRAY_DIM], dest[ARRAY_DIM][ARRAY_DIM];
  725. /*copy only filled sites, blank others, ie the blocked sites*/
  726. {
  727.     for (i = 0; i < (xnrow + 2); i++)
  728.         for (j = 0; j < (xncol + 2); j++)
  729.             if (source[i][j] == FILLED)
  730.                 dest[i][j] = FILLED;
  731.             else
  732.                 dest[i][j] = BLANK;
  733.     return;
  734. }
  735. /*END OF COPY_ARRAY            */
  736.  
  737.  
  738.  
  739. int work_test()
  740. /*move a pointer around the outside left edge of each sub-pattern*/
  741. /*starting at a filled site on the top row*/
  742. /*return true if reach bottom row, return false if no connection*/
  743. /*the arrays i_del and j_del hold a sequence of trial moves*/
  744. /*that maintains the pointer on the outside left edge of a pattern*/
  745. /*the pointer should either reach the bottom row or return to the top row*/
  746. /*if it returns to the top row, the search along the top edge continues*/
  747. /*mark trail of pointer with a string of TRAIL values*/
  748. {
  749.     for (j = 1; j < (xncol + 1); j++)
  750.     {
  751.         if (work_array[1][j] == BLANK || work_array[2][j] == BLANK)
  752.             continue;
  753.         k = 2;
  754.         i = 2;
  755.         work_array[1][j] = TRAIL;
  756.         work_array[2][j] = TRAIL;
  757.         while (i > 1)
  758.         {
  759.             if ((k_last = k % 4) == 0)
  760.                 k_last = 4;
  761.             for (k = (k_last - 1); k < (k_last + 3); k++)
  762.             {
  763.                 ii = i + del_i[k];
  764.                 jj = j + del_j[k];
  765.                 if (work_array[ii][jj] != BLANK)
  766.                 {
  767.                     i = ii;
  768.                     j = jj;
  769.                     work_array[i][j] = TRAIL;
  770.                     break;
  771.                 }
  772.             }
  773.             if (i == xnrow)
  774.                 return TRUE;
  775.         }
  776.     }
  777.     return FALSE;
  778. }
  779. /*END OF WORK_TEST            */
  780.  
  781.  
  782.  
  783.